package ua.lviv.franko.solvers.stressanalysis.oneDimention;

import java.util.ArrayList;
import java.util.Iterator;
import java.util.List;

import ua.lviv.franko.integration.ConstantFunction1D;
import ua.lviv.franko.integration.ConstantFunction2D;
import ua.lviv.franko.integration.FunctionWithoutDerivatives1D;
import ua.lviv.franko.integration.GaussQuadrature;
import ua.lviv.franko.integration.IFunction1D;
import ua.lviv.franko.integration.IFunction2D;
import ua.lviv.franko.integration.MatrixFunction1D;
import ua.lviv.franko.integration.Quadrature;
import ua.lviv.franko.integration.VectorFunction1D;
import ua.lviv.franko.matrix.RowType;
import ua.lviv.franko.matrix.hash.lu.HashMatrixLUStyle;
import ua.lviv.franko.matrix.hash.lu.LUFactorization;
import ua.lviv.franko.mlsFunctions.oneDimension.BaseFunction;
import ua.lviv.franko.mlsFunctions.oneDimension.Particle;
import ua.lviv.franko.mlsFunctions.oneDimension.weightFunctions.Gaussian;
import ua.lviv.franko.mlsFunctions.oneDimension.weightFunctions.IWeightFunction;
import ua.lviv.franko.visualisation.EFGResult1D;

public class TimoshenkoMembraneSimpleMFSolver {

	// general data
	protected static final int		QUADRATURE_POINTS_NUMBER	= 7;
	protected static final int		NUMBER_OF_UNKNOWNS			= 3;
	protected Quadrature			gau;
	protected EquationParameters	parameters;
	protected HashMatrixLUStyle		mainMatrix;
	protected List<Double>			bTmpVector;
	protected double[]				b;
	protected double				leftBoundery;
	protected double				rightBoundery;
	protected BoundaryCondition		leftBC;
	protected BoundaryCondition		rightBC;

	// problem data
	protected double				kKk1;
	protected double				kKk2;
	protected double				kKk3;
	protected IFunction1D			rho1;
	protected IFunction1D			rho3;
	protected IFunction1D			mu1;

	// data for base functions
	protected List<Particle>		particles;
	protected List<BaseFunction>	baseFunctions;

	public TimoshenkoMembraneSimpleMFSolver(final EquationParameters parameters, List<Particle> particles, double leftBoundery, double rightBoundery, int PolinomialDegree,
			IWeightFunction weight) {
		super();
		this.parameters = parameters;
		this.particles = particles;
		this.rightBoundery = rightBoundery;
		this.leftBoundery = leftBoundery;

		gau = new GaussQuadrature(QUADRATURE_POINTS_NUMBER);
		// gau = new TrapezoidalRuleQuadrature(QUADRATURE_POINTS_NUMBER);

		kKk1 = parameters.getE() * parameters.getH() / (1 - parameters.getV() * parameters.getV());
		kKk2 = parameters.getkPrime() * parameters.getG() * parameters.getH();
		kKk3 = parameters.getE() * Math.pow(parameters.getH(), 3) / (12 * (1 - parameters.getV() * parameters.getV()));
		rho1 = new FunctionWithoutDerivatives1D() {

			@Override
			public double calculate(final double xxx) {
				double h = parameters.getH();
				final double k1 = parameters.getK1().calculate(xxx);

				return gau.integrate(-h / 2, h / 2, new FunctionWithoutDerivatives1D() {

					@Override
					public double calculate(double x) {
						double f1 = parameters.getF1().calculate(xxx, x);
						return (1 + k1 * x) * f1;
					}
				});
			}
		};

		rho3 = new FunctionWithoutDerivatives1D() {

			@Override
			public double calculate(final double xxx) {
				double h = parameters.getH();
				final double k1 = parameters.getK1().calculate(xxx);

				return gau.integrate(-h / 2, h / 2, new FunctionWithoutDerivatives1D() {

					@Override
					public double calculate(double x) {
						double f1 = parameters.getF3().calculate(xxx, x);
						return (1 + k1 * x) * f1;
					}
				});
			}
		};

		mu1 = new FunctionWithoutDerivatives1D() {

			@Override
			public double calculate(final double xxx) {
				double h = parameters.getH();
				final double k1 = parameters.getK1().calculate(xxx);

				return gau.integrate(-h / 2, h / 2, new FunctionWithoutDerivatives1D() {

					@Override
					public double calculate(double x) {
						double f1 = parameters.getF1().calculate(xxx, x);
						return (1 + k1 * x) * x * f1;
					}
				});
			}
		};
		// creation of base functions;
		this.createBaseFunctions(PolinomialDegree, weight);
	}

	public void generateMatrix() {
		mainMatrix = new HashMatrixLUStyle(RowType.ARRAY_TYPE, (particles.size()) * NUMBER_OF_UNKNOWNS);
		// b = new double[(particles.size()) * NUMBER_OF_UNKNOWNS];
		bTmpVector = new ArrayList<Double>((particles.size()) * NUMBER_OF_UNKNOWNS);
		for (int i = 0; i < mainMatrix.getN(); i++) {
			bTmpVector.add(0.0);
		}
		for (int i = 0; i < baseFunctions.size(); i++) {
			final BaseFunction bfi = baseFunctions.get(i);
			for (int j = 0; j < baseFunctions.size(); j++) {
				final BaseFunction bfj = baseFunctions.get(j);
				MatrixFunction1D tmpFunc = new MatrixFunction1D() {

					@Override
					public double[][] calculate(double x) {
						double fiVal = bfi.calculate(x);
						double fiDer = bfi.calculateDerivaty(x);
						double fjVal = bfj.calculate(x);
						double fjDer = bfj.calculateDerivaty(x);
						double a1 = parameters.getA1().calculate(x);
						double k1 = parameters.getK1().calculate(x);
						return new double[][] {
								// first row
								{ kKk1 / a1 * fiDer * fjDer + a1 * k1 * k1 * kKk2 * fiVal * fjVal, kKk1 / a1 * k1 * fjVal * fiDer - k1 * kKk2 * fjDer * fiVal,
										-a1 * k1 * kKk2 * fiVal * fjVal },
								// second row
								{ -kKk2 * k1 * fjVal * fiDer + k1 * kKk1 * fjDer * fiVal, kKk2 / a1 * fjDer * fiDer + a1 * k1 * k1 * kKk1 * fiVal * fjVal, kKk2 * fjVal * fiDer },
								// third row
								{ -a1 * kKk2 * k1 * fiVal * fjVal, kKk2 * fjDer * fiVal, kKk3 / a1 * fjDer * fiDer + a1 * kKk2 * fiVal * fjVal } };
					}

					@Override
					public double[][] calculateDerivaty(double x) {

						return null;
					}

					@Override
					public int matrixSize() {
						return 3;
					}
				};

				double[][] tmp = gau.integrate(Math.max(Math.max(bfi.getXi().getXi() - bfi.getXi().getRho(), bfj.getXi().getXi() - bfj.getXi().getRho()), leftBoundery), bfi
						.getXi().getXi(), tmpFunc);

				double[][] tmp1 = gau.integrate(bfi.getXi().getXi(),
						Math.min(Math.min(bfi.getXi().getXi() + bfi.getXi().getRho(), bfj.getXi().getXi() + bfj.getXi().getRho()), rightBoundery), tmpFunc);

				for (int k = 0; k < tmp.length; k++) {
					for (int l = 0; l < tmp.length; l++) {
						// System.err.print(tmp[k][l] + "\t");
						tmp[k][l] += tmp1[k][l];
					}

				}

				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS, j * NUMBER_OF_UNKNOWNS, tmp[0][0]);
				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS, j * NUMBER_OF_UNKNOWNS + 1, tmp[0][1]);
				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS, j * NUMBER_OF_UNKNOWNS + 2, tmp[0][2]);

				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS + 1, j * NUMBER_OF_UNKNOWNS, tmp[1][0]);
				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS + 1, j * NUMBER_OF_UNKNOWNS + 1, tmp[1][1]);
				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS + 1, j * NUMBER_OF_UNKNOWNS + 2, tmp[1][2]);

				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS + 2, j * NUMBER_OF_UNKNOWNS, tmp[2][0]);
				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS + 2, j * NUMBER_OF_UNKNOWNS + 1, tmp[2][1]);
				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS + 2, j * NUMBER_OF_UNKNOWNS + 2, tmp[2][2]);

				// for (int k = 0; k < 3; k++) {
				// mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS + k, j *
				// NUMBER_OF_UNKNOWNS, tmp[NUMBER_OF_UNKNOWNS*k]);
				// mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS + k, j *
				// NUMBER_OF_UNKNOWNS + 1, tmp[NUMBER_OF_UNKNOWNS*k]);
				// mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS + k, j *
				// NUMBER_OF_UNKNOWNS + 2, tmp[NUMBER_OF_UNKNOWNS*k]);
				//
				// }
			}

			// right hand side
			double[] vector = gau.integrate(Math.max(bfi.getXi().getXi() - bfi.getXi().getRho(), leftBoundery),
					Math.min(bfi.getXi().getXi() + bfi.getXi().getRho(), rightBoundery), new VectorFunction1D() {

						@Override
						public double[] calculate(double x) {
							double fiVal = bfi.calculate(x);
							double h = parameters.getH() / 2.0;
							double k1 = parameters.getK1().calculate(x);

							double rho1Val = rho1.calculate(x);
							double rho3Val = rho3.calculate(x);
							double mu1Val = mu1.calculate(x);

							double sigma11Minus = parameters.getSigma11Minus().calculate(x);
							double sigma11Plus = parameters.getSigma11Plus().calculate(x);

							double sigma13Minus = parameters.getSigma13Minus().calculate(x);
							double sigma13Plus = parameters.getSigma13Plus().calculate(x);

							double sigma33Minus = parameters.getSigma33Minus().calculate(x);
							double sigma33Plus = parameters.getSigma33Plus().calculate(x);

							return new double[] { fiVal * ((1 + k1 * h) * sigma13Plus + (1 - k1 * h) * sigma13Minus + rho1Val),
									fiVal * ((1 + k1 * h) * sigma33Plus + (1 - k1 * h) * sigma33Minus + rho3Val),
									fiVal * (((1 + k1 * h) * sigma11Plus - (1 - k1 * h) * sigma11Minus) * h + mu1Val),

							};
						}

						@Override
						public double[] calculateDerivaty(double x) {

							return null;
						}

						@Override
						public int vectorSize() {

							return 3;
						}
					});
			bTmpVector.set(i * NUMBER_OF_UNKNOWNS, vector[0]);
			bTmpVector.set(i * NUMBER_OF_UNKNOWNS + 1, vector[1]);
			bTmpVector.set(i * NUMBER_OF_UNKNOWNS + 2, vector[2]);
		}

	}

	public BoundaryCondition getLeftBC() {
		return leftBC;
	}

	public BoundaryCondition getRightBC() {
		return rightBC;
	}

	public void setLeftBC(BoundaryCondition leftBC) {
		this.leftBC = leftBC;
	}

	public void setRightBC(BoundaryCondition rightBC) {
		this.rightBC = rightBC;
	}

	public EFGResult1D solve() {
		if (leftBC.isDirihlet()) {
			mainMatrix.addRows(3);
			BaseFunction bf = baseFunctions.get(0);
			double[] vals = bf.calculateAllFunctionsValuesOnly(leftBoundery);
			int size = mainMatrix.getN();
			for (int i = 0; i < vals.length; i++) {
				mainMatrix.setValue(size - 3, i * NUMBER_OF_UNKNOWNS, vals[i]);
				mainMatrix.setValue(size - 2, i * NUMBER_OF_UNKNOWNS + 1, vals[i]);
				mainMatrix.setValue(size - 1, i * NUMBER_OF_UNKNOWNS + 2, vals[i]);

				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS, size - 3, vals[i]);
				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS + 1, size - 2, vals[i]);
				mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS + 2, size - 1, vals[i]);
			}

			bTmpVector.add(leftBC.uVal);
			bTmpVector.add(leftBC.vVal);
			bTmpVector.add(leftBC.mVAl);
		}
		if (rightBC.isDirihlet()) {
			mainMatrix.addRows(3);
			BaseFunction bf = baseFunctions.get(baseFunctions.size() - 1);
			double[] vals = bf.calculateAllFunctionsValuesOnly(rightBoundery);
			int size = mainMatrix.getN();
			for (int i = 0; i < vals.length; i++) {
				mainMatrix.setValue(size - 3, (i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS, vals[i]);
				mainMatrix.setValue(size - 2, (i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS + 1, vals[i]);
				mainMatrix.setValue(size - 1, (i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS + 2, vals[i]);

				mainMatrix.setValue((i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS, size - 3, vals[i]);
				mainMatrix.setValue((i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS + 1, size - 2, vals[i]);
				mainMatrix.setValue((i + baseFunctions.size() - vals.length) * NUMBER_OF_UNKNOWNS + 2, size - 1, vals[i]);
			}
			bTmpVector.add(rightBC.uVal);
			bTmpVector.add(rightBC.vVal);
			bTmpVector.add(rightBC.mVAl);
		}
		b = new double[bTmpVector.size()];
		for (int i = 0; i < b.length; i++) {
			b[i] = bTmpVector.get(i);
		}
		// Matrix.printMatrix(mainMatrix);
		// double bigNum = 1000000000.0;
		// mainMatrix.setValue(0, 0, bigNum);
		// mainMatrix.setValue(1, 1, bigNum);
		// mainMatrix.setValue(2, 2, bigNum);
		// b[0] = leftBC.uVal*bigNum;
		// b[1] = leftBC.vVal*bigNum;
		// b[2] = leftBC.mVAl*bigNum;
		//
		// int size = mainMatrix.getN();
		// mainMatrix.setValue(size-3, size-3, bigNum);
		// mainMatrix.setValue(size-2, size-2, bigNum);
		// mainMatrix.setValue(size-1, size-1, bigNum);
		//
		// b[size-3] = rightBC.uVal*bigNum;
		// b[size-2] = rightBC.vVal*bigNum;
		// b[size-1] = rightBC.mVAl*bigNum;

		// for (int i = 0; i < mainMatrix.getN(); i++) {
		// for (int j = 0; j < mainMatrix.getN(); j++) {
		// System.out.print(mainMatrix.getValue(i, j) + "\t");
		// }
		// System.out.println();
		// // System.out.println(b[i]);
		// }
		// System.out.println("-------------");

		HashMatrixLUStyle.printMatrix(mainMatrix);
		LUFactorization lu = new LUFactorization(this.mainMatrix);
		double[] results = lu.calculate(b);

		// for (int i = 0; i < particles.size(); i++) {
		// System.out.println(results[i*NUMBER_OF_UNKNOWNS]);
		// }
		return new EFGResult1D(results, baseFunctions, particles);
	}

	private void createBaseFunctions(int polinomialDegree, IWeightFunction weight) {
		baseFunctions = new ArrayList<BaseFunction>();
		for (Iterator<Particle> iterator = particles.iterator(); iterator.hasNext();) {
			Particle particle = iterator.next();
			BaseFunction bf = new BaseFunction(polinomialDegree, particle, weight);
			bf.setNodeList(this.getParticleSupport(particle, weight));
			this.baseFunctions.add(bf);
		}

	}

	private List<Particle> getParticleSupport(Particle particle, IWeightFunction weight) {
		List<Particle> list = new ArrayList<Particle>();
		for (Iterator<Particle> iterator = particles.iterator(); iterator.hasNext();) {
			Particle p = (Particle) iterator.next();
			if (weight.intersectedSupport(particle, p)) {
				list.add(p);
			}
		}
		// return list;
		return particles;
	}

	public static void main(String[] args) {
		IFunction2D f1 = new ConstantFunction2D(1.0);
		IFunction2D f3 = new ConstantFunction2D(1.0);
		IFunction1D k1 = new ConstantFunction1D(0.0);
		IFunction1D a1 = new ConstantFunction1D(1.0);
		IFunction1D sigma11Plus = new ConstantFunction1D(0.0);
		IFunction1D sigma13Minus = new ConstantFunction1D(0.0);
		IFunction1D sigma11Minus = new ConstantFunction1D(0.0);
		IFunction1D sigma33Minus = new ConstantFunction1D(0.0);
		IFunction1D sigma33Plus = new ConstantFunction1D(0.0);
		IFunction1D sigma13Plus = new ConstantFunction1D(0.0);
		double v = 0.3;
		double e = 10000;
		double h = 0.5;

		EquationParameters eqP = new EquationParameters(a1, k1, e, v, h, f1, f3, sigma13Plus, sigma13Minus, sigma11Plus, sigma11Minus, sigma33Plus, sigma33Minus);
		double a = 0.0;
		double b = 1.0;
		int n = 10;
		double hStep = (b - a) / (double) n;
		double rho = 1.50000001 * hStep;
		List<Particle> particles = new ArrayList<Particle>();

		for (int i = 0; i <= n; i++) {
			particles.add(new Particle(a + i * hStep, rho));
		}

		TimoshenkoMembraneSimpleMFSolver s = new TimoshenkoMembraneSimpleMFSolver(eqP, particles, 0.0, 1.0, 1, new Gaussian());
		s.setLeftBC(new BoundaryCondition(0.0, 0.0, 0.0, true));
		s.setRightBC(new BoundaryCondition(0.0, 0.0, 0.0, true));
		// s.setRightBC(new BoundaryCondition(1.0, 1.0, 1.0, true));
		s.generateMatrix();
		// s.solve();
		s.solve().display(a, b);

	}

}
